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CHAPTER  5.  CONCLUSIONS 


CHAPTER  1 


INTRODUCTION 

Optimal  control  of  stochastic  systems  has  far  reach¬ 
ing  applications,  ranging  from  the  control  of  space  satel¬ 
lites  [1]  to  the  control  of  economic  systems  [2] .  It  is 
our  aim  to  show  general  results  applicable  to  many  differ¬ 
ent  systems.  Therefore  no  attempt  will  be  made  to  spe¬ 
cialize  our  results  to  any  one  specific  field. 

This  chapter  will  present  our  objectives,  as  well  as 
some  of  the  work  previously  done  in  this  area.  It  will 
also  present  an  overview  of  our  approach  and  a  brief  de¬ 
scription  of  the  system  model  used  throughout  the  sequel . 

1 . 1  Objectives 

This  work  proposes  an  approach  for  the  control  of  a 
general  class  of  discrete-time  noi. linear  stochastic  sys¬ 
tems.  The  system  model  incorporates  a  deterministic  lin¬ 
ear  portion  together  with  a  nonlinear  function  of  the 
state  and/or  control  vectors  in  combination  with  a  white 
noise  vector,  where  no  Gaussian  assumption  is  made. 

Under  certain  conditions  imposed  on  the  statistics  of 
the  additive  nonlinear  stochastic  term,  and  assuming  per¬ 
fect  state  information,  the  optimal  control,  which  mini- 
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mizes  a  quadratic  performance  index  subject  to  the  nonlin¬ 
ear  system  constraint,  is  shown  to  be  a  linear  function  of 
the  state  vector.  This  work  will  also  show  that  for  cer¬ 
tain  infinite  horizon  problems,  an  uncertainty  threshold 
can  be  found  such  that  the  designer  can,  a  priori,  put  an 
upper  bound  on  the  allowable  noise  covariance  to  obtain  a 
bounded  optimal  constant  feedback  control. 


1 . 2  Previous  Work 

The  overall  problem  of  optimal  control  of  dis¬ 
crete-time  nonlinear  stochastic  systems  is  quite  complex. 
Therefore,  much  previous  work  has  been  done  on  specialized 
cases . 

For  instance,  in  the  design  of  control  systems,  the 
exact  values  of  system  parameters  are  not  often  precisely 
known.  These  parameter  perturbations  can  be  caused  by 
environmental  effects,  operator  error,  equipment  aging 
etc.  Typically  these  parameter  variations  are  treated  as 
additive  stochastic  sequences,  but  in  some  cases  of  large 
variation  this  modeling  is  not  accurate.  In  these  in¬ 
stances,  one  is  left  with  modeling  these  uncertainties  as 
multiplicative  noise.  [3]  -  [8]  investigate  the  optimal 
control  and  stability  of  these  linear  discrete-time  sys¬ 
tems  with  multiplicative  noise.  These  systems  are  general- 
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ly  referred  to  as  bilinear  stochastic  systems,  systems 
with  multiplicative  noise,  or  stochastic  parameter  sys¬ 
tems.  An  excellent  review  of  the  continuous-time  counter¬ 
part  of  this  class  of  systems  is  given  in  [9] . 

[10]  -  [12]  deal  with  the  estimation  problems  related 
to  stochastic  parameter  systems.  These  problems  are  asso¬ 
ciated  with  phenomena  such  as  fading  or  reflection  of  a 
transmitted  signal  from  the  ionosphere,  and  certain  situa¬ 
tions  involving  sampling,  gating  or  amplitude  modulation. 

The  more  general  case  of  discrete-time  nonlinear 
stochastic  systems  is  investigated  in  [13] ,  which  was  used 
as  a  basis  for  much  of  our  work.  [14]  and  [15]  suggest 
several  novel  control  schemas  for  these  systems. 

[16]  -  [18]  address  the  steady  state  characteristics 
of  linear  stochastic  parameter  systems.  [16]  derives  sta¬ 
bility  conditions  for  the  matrix  Riccati  equation  arising 
in  the  optimal  control  of  linear  systems  with  random  gain. 
[17]  and  [18]  derive  an  uncertainty  threshold  for  the 
existence  of  the  infinite-horizon  solution  to  the  optimal 
control  of  linear  discrete-time  stochastic  parameter  sys¬ 
tems.  This  work  is  generalized  to  nonlinear  stochastic 
systems  in  chapters  3  and  4 . 
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1 . 3  Overview  of  Our  Approach 

Our  intent  is  to  design  the  optimal  linear  feedback 
controller  for  a  general  stochastic  system  configuration. 
We  will  endeavor  to  keep  the  system  as  general  as  possible 
so  that  the  resulting  controller  will  be  suited  for  as 
wide  a  range  of  applications  as  possible.  In  the  remain¬ 
der  of  this  chapter  we  describe  this  general  system  and 
then  give  several  examples  of  both  linear  and  nonlinear 
systems  that  it  covers.  The  associated  optimal  fi¬ 
nite-horizon  controller  for  this  system  will  be  presented 
in  chapter  2. 

The  general  finite-horizon  control  solution,  although 
quite  interesting  academically,  is  not  always  the  most 
practical  approach.  Quite  often  it  is  advantageous  to 
utilize  a  constant  feedback  controller.  This  prompts  us 
to  investigate  the  infinite-horizon  solution.  Our  concern 
here  is  does  the  infinite-horizon  solution  exist?  And,  is 
there  a  quantifiable  means  to  predetermine  this  existence? 
Through  steady-state  analysis  of  a  Riccati-like  equation 
in  Chapters  3  and  4,  threshold  conditions  are  developed 
for  both  the  scalar  and  the  multivariable  cases,  such  that 
the  designer  can,  a  priori,  based  on  the  covariance  of  the 
stochastic  parameters,  guarantee  the  existence  of  a 
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steady-state  constant  feedback  controller.  Computer  simu¬ 
lations  follow  both  developments  to  verify  these  threshold 
conditions . 


1 . 4  General  System 


The  optimal  control  problem  is  to  minimize  the  fol¬ 
lowing  quadratic  performance  index 

e{£xT"s*x"+  +  +  (1.1) 


where  N  is  the  final  time,  ZT(  .  >  denotes  the  expectation 
operator,  S*  and  Qk  are  positive  semi-definite  symmetric 
matrices,  Rk  is  a  positive  definite  symmetric  matrix,  and 


ft 


subject  to  the  system  constraint 


+  Bkut  +  fk[xk,uk,at) 


(1.2) 


with  x0  given,  and 


xkeR*,  ukeR",  a teR*,  fk:RMxRmxR*-*R * 


The  noise  sequence  a,  is  assumed  to  be  independently  dis¬ 
tributed  in  time,  but  not  necessarily  Gaussian,  and 
/»{xk,u»,at)  has  the  following  statistical  description 


E {/*(**»u*»a*)l**u’*}“0  V  xkeR*,  ukeRm,  fc-0 N-\ 


(1.3) 


f  *(*»»“*)*  £{/*(**. «»)/!(**. |  xtu,}  (1.4) 

is  a  quadratic  function  of  xt  and  uk  having  the  following 


form 


F t[xk,iik)^  Pl+^P^xlW'tXt* xlN'ku.k  +  }^u.lhilku^  (1 .5) 


where 


Pkt  W't-.  and  Mk  are  symmetric  and  n'  ■  ^ 

J+ 


Fk[xktuk)ZO  V  x.Eft*, 

Note  that  the  above  condition  is  not  at  all  restrictive  as 
it  is  necessary  in  order  for  equation  (1.5)  to  be  a  proper 
covariance  representation. 


i  r*  ri-i 


r*.T 


v\nV 
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Description 


The  above  general  problem  description  encompasses 
many  known  standard  systems.  For  instance,  one  of  the 
most  common  examples  is  a  linear  system  with  an  additive 
noise  term.  Let 


ft(xkt  ut,aj-ak 


thus  reducing  equations  (1.2)  -  (1.5)  to  the  following 


**.i  -  AtXi  +  Bku.t+  a, 


(1.6) 


with 


(1.7) 


Similarly  (1.2)  -  (1.5)  can  be  used  to  represent 

linear  stochastic  parameter  systems.  The  general  form  is 


-  At(<v)xk  +  Bt(u>)ut 


(1.8) 


where  xk(a>)  and  5»(w)  are  matrices  having  elements  which  are 
white  noises  possibly  correlated  with  each  other  at  each 
time  instant  k .  This  modeling  evolves  naturally  in  sam¬ 
pled  versions  of  diffusion  processes  associated  with  nu¬ 
clear  fission  and  heat  transfer,  as  well  as  in  the  migra¬ 
tion  of  population,  and  the  growth  of  biological  cells, 


etc.  This  is  also  the  model  typically  used  to  describe 
the  uncertainties  in  the  parameters  connected  with  econom¬ 
ic  predictions  [2].  One  can  also  arrive  at  Cl. 8)  by: 


1.  Uniformly  sampling  a  continuous-time  stochastic 
parameter  system. 

2.  Randomly  sampling  a  continuous-time  deterministic 
syst em . 

3.  Modeling  a  deterministic  system  with  an  additive 
noise  having  random  correlation  characteristics. 

The  first  two  are  self  explanatory  and  the  third  is  shown 
below.  Suppose 


xk*i  m  Akxk  + Bkuk  + Ckak  (1*9) 

with  the  noise  term  at  possessing  the  dynamics 

at.l-Dt(aj)ak  +  0k  (1.10) 

where  0»(<u)  is  a  matrix  with  random  elements,  and  fit  is 
additive  white  noise.  Then  by  enlarging  the  state  space, 
we  obtain 


■**♦1  Ak  Ck  xk  +  Bk  0 

"L°  0*(<u)JLa*J  L  o  j  *  +  L£*_ 


(l.ii) 


iV*'*1****  i*«  •*»  **. 
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1 . 6  Example*  of  Nonlinear  Stochastic  Systems 


The  following  are  some  novel  special  cases  of  the 
general  system  described  by  equations  (1.2)  -  (1.5), 


1 .  Norm  Dependent  Random  Vector  : 


Then 


with 


/*(**.  u*,aJ-at^Af[Z>i*k  +  x[D2u*  +  ^uf£>3uk  (1.12) 


.  £{«*}-<> 


O3>0  and  [D,  -  D2D'9'Dt2]*  0 


(1.13) 


f*(**.ttk)-/>k^xlD,xk-»-xlD2ttk  +  |a[D3UkJ  (1.14) 


Note  that  if  D{- 1,  i>j  =  ,Di  =  0,  then 


/*(«». ak,ak)-ak||xk 


(1.15) 


In  the  following  examples  only  the  system  will  be 
given,  as  the  covariance  matrix  can  be  determined  by  fol¬ 


lowing  the  above  procedure. 


2.  Random  Vector  Dependent  Upon  the  Sign  of  a  Scalar 
Nonlinear  Function  of  xk  and  uk : 

Let 


/«(**. u*,at)-  sgn  [<#>(**  ,uk)]at  (1.16) 

where  ,  and  the  statistics  of  ak  are  as  in  equa¬ 

tion  (1 . 13)  . 

3.  Random  Vector  Dependent  Upon  the  Absolute  Value  of  a 
Linear  Combination  of  xk  and  uk 

Let 


/*(**.  +  (1.17) 

Again  the  statistics  of  «k  are  given  in  (1.13) . 

4 .  Random  Vector  Dependent  Upon  the  Norm  of  xk  and  the 
Absolute  Value  of  Components  of  xk 

Let 


(1.18) 


where  the  n  vectors  . <*;  are  uncorrelated  and  »l  is  the 

t"  component  of  xk  . 
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It  is  worth  noting  that  the  above  examples  are 
uinely  nonlinear,  and  no  assumption  has  been  made  on 
type  of  probability  distribution  of  at  . 
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CHAPTER  2 


OPTIMAL  CONTROL 


This  chapter  presents  the  optimal  solution  to  the 
general  discrete-time  nonlinear  stochastic  system  outlined 
in  Chapter  1  equations  (1.1)  -  (1.5)  .  Our  approach  is  to 
merely  state  the  solution  now,  and  then  in  section  2.4  use 
stochastic  dynamic  programming  to  prove  that  this  solution 
is  in  fact  optimal.  Rather  than  stating  the  solution  di¬ 
rectly,  we  will  first  try  to  gain  insight  by  examining 
known  simpler  cases. 


2 . 1  Deterministic  Case 


The  deterministic  case  is  well  documented  in  most 
texts  on  optimal  control.  [19]  is  one  such  reference. 
Once  again  the  performance  index  is  a  quadratic  function, 
but  there  is  no  need  to  take  the  expected  value,  due  to 
the  deterministic  nature  of  the  problem 


1  j  w-i 

~  N  +  -  £  (x l Q kx  *  x l *  aTt Rtut ) 

*•  *»o 


(2.1) 


where  sw  and  are  positive  semi-definite  symmetric  ma¬ 

trices,  Rt  is  a  positive  definite  symmetric  matrix,  and 
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The  system  model  is 


**♦1-  AtXk  +  BkUi 


(2.2) 


For  this  system  the  optimal  control  is  known  to  be 


(Xu  m  ~K  kXk 


(2.3) 


where  Kk  is  called  the  Kalman  gain,  which  is  given  by 


(2.4) 


J»  is  the  solution  backward  in  time,  from  5„ ,  of  the 
following  Riccati  equation 


(2.S) 


And  the  optimal  cost  is  given  by 


J |  "  Stx ( 


2.2  Stochastic  Parameter  CaBe 


(2.6) 


References  [4]  and  £18]  make  use  of  the  following 
system  description.  We  find  this  case  interesting  because 
it  is  very  similar  to  the  general  form  to  be  proposed  in 


k«JLU, 


section  2.3.  Also,  much  of  the  derivation,  although  not 
presented  here,  is  the  same.  We  minimize  the  following 
performance  index 


e{^xtmS Mx Y.  ^^xlQtXk  +  xlTtUt+  (2.7) 

where  S0  and  Q*  are  positive  semi-definite  symmetric  ma¬ 
trices,  A,  is  a  positive  definite  symmetric  matrix,  and 


subject  to  the  following  stochastic  system  model 


•**♦1  "  <*kAkXk  +  6kBkU.k 

(2.8) 

where  at 

with  the 

and  a,  Are  scalar  Gaussian  white  random 

following  known  stationary  statistics 

sequences 

£(a,|-5.  £{(o,-5)2}-r 

(2.9) 

E(«,)-5,  £{(«,-5)2)-/) 

(2.10) 

(2.11) 

The  form 

of  the  optimal  control  is  once  again 

Ukm-KkXk 

(2.12) 

But  K t  is  now  given  by 


/rt-[(52  +  j)fij[S)k.15lt  +  Rt]‘,[Ca5+/l)flISt*1^»  +  r[]  (2.13) 

and  s»  is  given  by 

S t  m  (a2+  r)AlSt.t/4 1  +  Qt  + 

~[(a8  +  A)AZSt.lBlt  +  Tlt}[(82  +  A)BlSk.lBk  +  Jlt]',l(a5+A)BiSt.lAk  +  Tl] 

(2.14) 

Notice  the  similarity  between  equations  (2.4)  and  (2.13) 
and  between  equations  (2.5)  and  (2.14).  Also  the  expres¬ 
sion  for  the  optimal  cost  remains  unchanged. 


(2.15) 


J ,  *  —x  i  S fX  t 


2 . 3  General  Nonlinear  Stochastic  Case 


The  optimal  solution  to  the  system  descxibed  in  sec¬ 
tion  1.4  is  as  follows.  The  optimal  control  has  the  same 
feedback  form 


ut« -Kkx, 


(2.16) 


but  AT,  is  now  given  by 


*.-[*.♦  BTtSk.,Bt*  Jr.]"(T['.Si[S..1/1..Wr]  (2.17) 


St  is  given  by 


S  tm  Qt*  AlSk.tAk  +  -  +  W  k]K  k  (2.18) 
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and  the  optimal  cost  is  given  by 


J  im^x\S  tXt*Qt 


(2.19) 


where  5»  is  given,  ««s0  and 


a*  "  2  8‘*> 


(2.20) 


^ -  iSil  tr 


(2.21) 


W,-rII  tr  (s,., />;)*;] 

^<•1 


(2.22) 


(2.23) 


2.4  Proof  of  the  General  Result 


The  proof  of  the  general  result  can  be  broken  down 


into  three  main  steps.  First,  stochastic  dynamic  program¬ 


ming  or  Bellman's  principle  of  optimality  [19]  will  be 


introduced.  Then,  using  this  concept,  we  will  generate  a 


form  of  the  Bellman  functional  equation  for  the  optimal 


cost  to  go  from  time  k  to  finish  time  N.  Finally,  induc¬ 


tion  will  be  used  to  verify  a  proposed  solution  to  this 


-16- 


functional  aquation.  The  solution  for  the  general  nonlin¬ 
ear  stochastic  case,  equations  (2.16)  through  (2.23),  is  a 
result  of  this  induction  step. 

i)  Stochastic  Dynamic  Programming 

With  dynamic  programming,  we  are  not  concerned  with 
how  we  came  to  be  in  the  present  state,  but  are  only 
concerned  with  optimizing  all  future  decisions,  which  is 
called  Bellman's  principle  of  optimality  [19].  If  this 
procedure  were  carried  out  iteratively,  the  result  would 
be  an  overall  optimal  control  sequence.  With  this  thought 
in  mind,  let  us  examine  equation  (1.1),  which  is  repeated 
here  for  convenience 

(i.l) 


Equation  (1.1)  may  be  written  in  two  parts  (notice 
the  temporary  change  in  subscripts) . 


Q tX l  +  x  *  T t^’t*  2U"i 

■  E  ly~x  MS  tjX  M  +  ^^2*  t  Q  *+ x  t  T  tU-t* —u,(  (2.24) 


We  are  trying  to  minimize  equation  (2.24)  by  appropriately 
choosing  the  control  sequence  u»  .  The  first  term  of  equa- 
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tion  (2.24)  does  not  depend  on  uk  ,  ,  <*».»,  .  ..  >  “*-•  • 

Therefore,  applying  Bellman's  principle  of  optimality,  it 
is  sufficient  to  minimize  only  the  second  term  over  u»  , 

***♦■  >  “*♦*  >  •  •  •  >  u»-i  • 


ii)  The  Bellman  Equation 


In  analyzing  the  second  part  of  equation  (2.24),  the 
following  properties  prove  useful  and  can  be  found  in 
[21].  Here,  Y  is  an  arbitrary  function  of  x. 


E{Y}-E{E[Y\x]} 


where  the  right  hand  side  of  the  above  equation  is  the 
expected  value  of  the  conditional  expectation  of  Y  given 
x.  Also 


min£<y}-£{miny> 


putting  these  together 


min  £{/>  -  E  {min  E[Y  |  x]> 


Applying  all  this  to  the  right  hand  side  of  equation 
(2.24)  we  have 

min£ ^xlS„xM+  j  {^xjQtx,+  xjT tut+  -  £{j  »(**)} 

(2.25) 
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where 


Jk(xj- jnin  e1^xt„Shxh  +  £  (^*fQ<*«  +  *f7\u<+iufrttu<j|x*j 


(2.26 


Representing  J*(xk|  as  a  summation  of  minimizations  over  u» 
U|«i  }  U|4t  f  •  •  •  9  tip- 1  f  we  Have 


J  k(xk)~minE!^xlQkxk  +  xlTtut  +  ^ulRkut- 


+  m 


*♦1 


+  min  —  Xt*^Qk*2Xk»2*  Xk*2^"  k»2Uk*2*  2a*»2^  k*2U  k»3 


min  E...  j  x,.2 


jl**j 


(2.27 


One  can  easily  see  that,  excluding  the  first  term,  aqua 
tion  (2.27)  is  simply  by  our  definition  of  J,(x») 
Substituting  J  )  into  (2.27)  results  in 


A(x*)  -  nrun ^E^xrtQkxk*  xTtT kut+  ^u£rttut  |  x,J  +  £p*.,(***i)  I  **} ] 


(2.28 


As  u*  is  a  function  of  xt ,  all  the  quantities  in  the  first 
expected  value  of  equation  (2.28)  are  known  constants,  for 
a  given  ft.  This  results  in  the  following  form  of  the 
Bellman  functional  equation 

•/ *(**)-  +  I  **}] 


(2.29) 


iii)  Solution  of  the  Functional  Equation 


The  solution  to  equation  (2.29)  will  be  shown  to  be 


A(*J“  +  S*  given,  e„m0 


(2.30) 


First,  by  definition 


J n[x 2 *TnS nx N 


(2.31) 


Next ,  assume 


J  t*  i  ( x  »♦!  )  “  q  x  ** 1  $  ** 1  x  *♦ 1  +  a  1 


(2.32) 


Finally,  to  show  that  equation  (2.32)  implies  equation 
(2.30),  we  will  start  by  inserting  equation  (2.32)  into 
equation  (2.29) 
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J*(x J- + xlT^u*  +  iujA,u*+ +  «*.,  |x*jj 


(2.33) 


Inserting  equation  (1.2)  for  x»., 


Jt{xt)mmin^xlQkxk*xlTkak*^alRkat* 

+  f|— [yitXlt+ fifcU,k  + /k]  5lk*l[^tXfc+  BkUk  +  /*]+  ®fc*i 


(2.34) 


Expanding  equation  (2.34) 


«/*(*«)“ mjn^J£[<3»JCk  +  x[Tkuk  +  ^u[Rltult  + 

+  t*iAkxk  +  —U^B^S  n.xBiUi*  —fkS  mf  t 


+  — Xfc>l45fc,iflJku(t+— u»5jSn*|>1fcxt  +  ~Xnilj5k,i/t  + 


+  2^*  ^  *  $  *41  +  ^  /t*  *  +  f  *  ^**i  ^  *u*  +  ®**i  l**|^  (2.35) 


Now,  let  us  consider  equation  (2.35)  in  depth.  First,  the 
expected  value  can  be  distributed  over  the  sum  of  terms . 
Second,  because  x»  is  known  (assuming  complete  state  in- 


pu  vl  Arwwvim'www¥uw^n#j'ir.'  /  v.vjwwm  tvi  a*  rum  ah  ru*  r  j<  n*  vn^n/v  'UHrjr^y  wjr^^jrF?r*jr»^r’-v’-jr\KT. 


formation),  all  of  the  terms,  with  the  exception  of  terms 
containing  /, ,  can  be  considered  constants  for  a  given  k . 
Furthermore,  the  expected  value  of  any  term  containing  a 
single  /»  term  is  zero,  as  the  expected  value  of  /*  is 
defined  to  be  zero.  Also,  each  term  is  a  scalar  (lxl). 
Therefore,  each  term  is  equal  to  its  transpose,  allowing 
us  to  further  reduce  the  equation 

^*(**)“  min  [  l;xlQkxk  +  xrkT  ku.k  +  +  + 

+  **i & kat +  u’l&lSk*iAkxk+  B^—/lSk^.1/k  I  JC*| +  j 

(2.36) 


Now  we  need  to  examine  the  remaining  expected  value  term 
in  equation  (2.36).  As  was  stated  earlier,  we  are  working 
with  a  scalar  term.  This  allows  us  to  take  the  trace  of 
the  inside  term,  as  the  trace  of  a  scalar  is  equal  to 
itself.  Also,  the  trace  and  expected  value  may  be  inter¬ 
changed,  and  the  expected  value  of  $*,  =  $*.,  .  Using  these 
facts  and  the  property 


tr  [ABC]  -  tr  [CAB] 


(2.37) 


yields  the  following  equations 


WWW 


I **}-*{!  tr  tr  [/*/j^I]l*,}- 

“^2  I  x*| "  2  *r  I  x*l] “  2  *r  I  *»}] 

(2.38) 

Looking  at  the  rightmost  equation  of  (2 . 38) ,  one  can  see 
that  the  expected  value  is  simply  Ft  from  equations  (1.4) 
and  (1.5).  Therefore,  substituting  equations  (1.5),  and 

(2.38)  into  equation  (2.36)  yields 

+  2  u’*BtSt.iBtu.it  +  u.%B%St.lAkxk  +  6k.l  + 

tr  tr  + 

+  5  tr  tr  [5**‘Z/>*U*V*“*J 


(2.39) 


We  need  to  investigate  the  last  three  terms  further,  and, 
as  they  are  basically  the  same,  we  will  consider  just  the 
first  of  the  three.  The  trace  and  can  be  taken  inside 
the  summation,  as  shown  in  equation  (2.40). 

t  tr  I  Sw.f  tr  (2.40) 


VmVViVi 


5*m*o 


Using  the  fact  that  x;v‘.x*  is  a  scalar,  we  can  remove  it 
from  the  trace  leaving 


(2.41) 


Now  we  can  insert  the  trace  between  *»  and  w, ,  as  the 
trace  is  simply  another  scalar.  Also  xk  does  not  depend 
on  i  and  can  be  removed  from  the  summation. 


1 


r 

t 


S 


t*i 


^xTkWkxt 


(2.42) 


where 


17, 


(2.43) 


Similar  procedures  can  be  followed  to  simplify  the  last 
two  terms  in  equation  (2.39),  generating  the  following  two 
definitions 


tr  [S,.,Pi)wL  (2.44) 

x(ml 

tr  (2.45) 

l 


Substituting  all  of  this  into  equation  (2.39),  and  group¬ 
ing  like  terms  results  in  the  following  equation 
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*/*(«*)- mi n  1^ ,+ /l  JS + 17 (tjxk+xrrJka/t  +  -arRfcaJk  + 


+  2  B  t  S  i*i  B  ku.t  *  u\B  kS  fc»i  A  kx  t  +  fl*»i  + 


2  ^"r  *1 +  x * +  2 ^ j  (2.46) 


To  find  the  minimum,  we  can  take  the  partial  derivative  of 
equation  (2.46)  with  respect  to  u»  ,  and  set  this  equal  to 


J *lx t) m  T%x k  •*•  Rkuk‘*‘  B kS k4.t  B kuk  +  B%S k.t  A kx k  +  N  t  x  t  +  W  ku.k  m  0 

(2.47) 


grouping  like  terms 


lRk  +  BTkSk.lBk+IIk]uk--[TTkxk  +  BTkSk.lAkxk  +  nTkxk]  (2.48) 

Vt  and  Sf i  can  be  shown  to  be  positive  semi-definite  and 
positive  definite  respectively,  and,  since  Jt»  is  positive 
definite  the  following  inverse  exists 


[Rk  +  BlS^Bt  +  IT,]-1 


Therefore,  solving  for  u» 


u(t--[A^fl[StM5t+W't]-,[Tl  +  aIS4M^t+77l]x,  (2.49) 
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As  can  be  seen  in  equation  (2.49),  we  have  finally  arrived 
at  the  optimal  controller  proposed  in  equations  (2.16)  and 
(2.17).  We  have  also  shown  the  derivation  of  equations 
(2.21)  -  (2.23),  but  it  remains  yet  to  complete  the  induc¬ 
tion  on  Jt  ,  and  verify  equations  (2.18)  -  (2.20).  This 

can  be  done  by  substituting  equation  (2.49)  into  equation 
(2.46).  After  appropriate  algebraic  manipulations  this 
results  in 

J  k(xk) m  +  ^  k  ~  *]^  *}•**■*■ 

+  \  tr  (2.50) 

We  began  the  induction  step  with  J and  needed  to  show 
that  this  implied 


J *(*•)“  ^*TkSk*k  +  *k 


(2.51) 


Comparing  equation  (2,51)  with  equation  (2.50)  we  can  see 
that  this  is  true  if  we  assume  the  following 


Sk-Qt+ATlSk.lAt  +  W't-KHRt  +  BlSk.lBk+ITt]Kt  (2.52) 


6kmn  + 


(2.53) 
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r  +  St*,(ab  + j) 

. - 7 - — r  x  * 

ft  +  S**,(b2+7) 


u.  -  - 


(3.4) 


.  ,0,5  f a,.V\  [T»5,.l(ab  +  f  )]* 


(3.S) 


with  the  optimal  cost  to  go  from  time  I  equal  to 


(3.6) 


3.2  Threshold  Condition 


Consider  equation  (3.5)  and  suppose  5^,  gets  very 
large.  Then  we  can  consider  Q,  T,  and  R,  added  to  a  very 
large  $*.,  ,  as  negligible,  thus  reducing  equation  (3.5)  to 


*  **‘V 


(3.7) 


As  we  are  only  interested  in  the  steady-state  value  of  S»  , 
forward  or  backward  evolution  in  time  does  not  matter, 
therefore,  we  will  interchange  5>  and  s*.i . 


S*.x  " 


(ab 


5. 


(3.8) 


or 
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5t*(  *  g  5fc 


where 


•  -(-a-fep 


(3.9) 


Equation  (3.9)  is  the  threshold  condition.  It  is 
obvious  that  if  g  <1  then  the  evolution  of  Sk  stays 
bounded  and  S*  ,  and  subsequently  K, ,  has  a  steady-state 
solution . 


There  is  an  alternate  approach  to  finding  the  thresh¬ 
old  condition  g  that  does  not  utilize  the  above  approxima¬ 
tions.  First  suppose  st  does  reach  a  steady-state,  say  5, 
then  equation  (3.5)  takes  on  a  somewhat  simpler  form. 


„_c.  J\a,  [r*j(«t*5)r 

S  5  Sl°  2)  -**s(6**S) 


(3.10) 


multiplying  both  sides  by  the  denominator  and  rearranging 
terms  yields 


[7*5(o6*f)]  -[ft*5(6s*|)][Q*s(aJ*!£)-s] 


(3.11) 


Next,  squaring  the  left  side  and  multiplying  out  the  right 
side  yields 
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T2  +  2TS^a6  + jj  +  52^ab  + jj2-RQ-RS^a2+ jj  +  RS- 

-«»(fc**|)*S*(fc**f)(«**K)-S*(6**|)  (3.12) 


finally  combining  the  like  powers  of  5 


■5[/?(<X2+j)‘27'(ab  + j)  +  <?(ba+f)“A]  +  /?0"r2“0 


(3.13) 


Upon  examining  equation  (3.13)  one  can  see  that  it  is  just 
a  quadratic  equation  in  S,  of  the  general  form  shown 
below. 


AS2*  BS  +  C  -  0 


(3.14) 


Thinking  in  two  dimensions,  if  A<o,  then  we  have  a 
parabola  concave  downward.  If  at  the  same  time  c>o,  then 
the  Y  intercept  is  positive  and  we  have  a  unique  positive 
definite  solution  for  s.  Now,  let's  investigate  these 
conditions  with  respect  to  equation  (3.13).  C> o  implies 
*Q-r*»o.  This  is  a  natural  restriction  as  T  is  the  cross 


weighting  term  from  x»  and  ii*  and  can  not  be  weighted 


larger  than  Q  or  jt,  which  are  the  weighting  terms  associ¬ 
ated  with  *;  and  respectively.  As  a  result  AQ-r*>o  is 

automatically  true. 

Therefore  all  that  is  necessary  for  5  to  have  a 
unique  positive  definite  solution  is  the  following 

[(b2t0“2*T)-(a6tl)2-(!,2*f)]<o  (3'15) 

Simplifying  equation  (3.15)  we  arrive  at  the  threshold 
condition 


V  2)  b2  +  Z 


(3.16) 


It  is  worth  noting  that,  as  the  above  derivations 
show,  the  threshold  condition  g  is  both  a  necessary  and 
sufficient  condition.  If  g  is  less  than  one,  then  the 
steady-state  solution  exists,  but,  if  g  is  greater  than 
one,  then  5k  increases  without  bound. 
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3.3  Simulations 

^  1 
•'V 

Figures  3.1 

through  3,3  plot 

the  evolution 

of  5  k 

y 

given  by  equation 

(3.5) ,  forward  in 

time,  vs  k  for 

fifty 

■\ 

► 

values  of  k .  Each 

graph  contains 

four  plots  of 

5k  for 

I 

various  values  of 

V  ,N,  and  u  .  These  values  were 

specif- 

<  •  •’ 
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ically  chosen  to  demonstrate  the  threshold  condition.  It 
should  be  noted  that,  as  the  threshold  condition  is  ap¬ 
proached,  the  values  of  5«  become  very  large.  Therefore 
the  base  ten  logarithm  of  S„  vs  t  is  what  actually  appears 
in  the  graphs . 

It  should  also  be  noted  that,  as  we  are  plotting  St 
forward  in  time,  the  equation  for  the  optimal  cost  to  go 
from  time  0  to  time  k  is  now 

Sk  (3 .17) 

As  shown  in  equation  (3.17),  the  optimal  cost,  for  a  given 
x«  ,  is  merely  a  constant  times  S*  .  Therefore,  figures  3.1 


through  3.3  can  also  be  considered  the  evolution  of  the 


Al 

f 


EVOLUTION  OF  Sk 

SCALAR  CASE 


wwhjw  rw  njf  w>  v*  wwwwuvw  \^7V;  '/vi-'w. 1 


CHAPTER  4 

E  INFINITE-HORIZON  CASE 


Unlike  the  scalar  threshold  condition,  equation 
(3.16),  an  absolute  threshold  for  the  multivariable  case 
is  much  more  difficult  to  define.  Therefore,  a  sufficient 
condition  is  found  that  guarantees  the  existence  of  the 
infinite-horizon  solution.  If  the  threshold  is  violated, 
the  steady-state  solution  may  or  may  not  exist.  In  spite 
of  this,  the  relative  magnitude  that  the  threshold  condi¬ 
tion  is  exceeded  by  is  still  an  excellent  indicator  of  the 
steady-state  characteristics  of  Sk . 


As  before,  we  are  only  interested  in  the  steady-state 
behavior  of  5k ,  and  we  can  choose  to  evaluate  this  forward 
in  time  rather  than  backward.  Therefore  $k ,  which  is 
equation  (2.18),  will  be  rewritten  forward  in  time.  As  we 
are  interested  in  the  steady-state  analysis,  all  matrices 
except  Si,  will  be  assumed  constant.  Also,  to  make  the 
problem  more  tractable,  without  a  great  loss  of  generali¬ 
ty,  we  will  drop  the  cross  terms  rk  and  W*  and  assume  Q  is 
positive  definite.  Therefore  we  will  begin  with  the  fol¬ 


lowing  equation 


I  *- 


n 


Sk.t-Q  +  ATSkA  +  ^Y.  tr  (SS*)*'* 

<*  1-1 


->lr5<5^R  +  BrS*fl  +  ^|j  tr  BT S kA  (4.4) 


We  will  now  show  three  properties  of  the  evolution  of 
5k  that  guarantee  the  existence  of  a  steady-state  solu¬ 
tion.  First,  we  will  show  that  St  is  always  positive 
definite.  Next  we  will  show  that  St  is  a  monotonically 
increasing  matrix  sequence.  And  finally,  if  we  also  show 
that  5t  remains  bounded  above,  then  it  must  come  to  a 
steady-state  value  [22] .  Therefore,  we  will  prove  the 
first  two  properties,  and  then  find  a  threshold  condition 
such  that  the  third  is  true. 
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First,  we  will  show  that,  starting  from  s»- 0,  all 
subsequent  5*  are  positive  definite.  Evaluating  equation 

(4.4)  with  *=0  yields  and,  as  <3  is  positive  defi¬ 
nite,  Si  must  be  positive  definite.  Now  we  may  apply  the 
well  known  matrix  inversion  lemma  [19] 


[vl  II  +  >1 12 >1  22-^21  ]  m  ^  11  “  A  11  A  t2[^21  A  n  A  u  +  ^22]  •^21^11 


where  A,t  and  An  are  square  and  non-singular,  to  equation 

(4.4)  to  obtain 

tr  arj  A+Q+'-f'tT  [SkPl)Wl 

(4.5) 

Substituting  S,=Q,  which  is  positive  definite,  into 
equation  (4.5),  one  can  see  that  the  two  summation  terms 
are  at  least  positive  semi-definite,  and,  as  R  is  positive 
definite,  the  inner  inverse  exits  and  is  positive  defi¬ 
nite.  Now,  using  these  facts  and  the  property  that  any 
matrix  quadratic  form  zxzT  is  at  least  positive  semi-defi¬ 
nite  if  X  is  at  least  positive  semi-definite,  we  can  show 
that 

+  tr  BT*0 
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We  have  shown  s,  is  positive  definite,  therefore  s;'  ex¬ 
ists,  and,  moreover,  s;‘  is  positive  definite  which  means 
the  outer  inverse  exists  and 


+  +  tr  (SfcP‘)A/‘j  A  *  0 


We  already  have  <2  positive  definite.  Therefore,  St  is 
positive  definite.  Obviously  this  pattern  continues  with 
increasing  Ik.  Therefore,  starting  from  $.=0,  all  subse¬ 
quent  St  are  positive  definite. 


Next  we  will  show  that  Sk  is  a  monotonically  increas¬ 
ing  matrix  sequence.  We  have  already  shown  that  starting 
from  5«=0,  S,=Q,  which  is  positive  definite.  Using  this 
as  the  basis,  we  will  now  show  that  5k25».|  implies  Jk.,i^k  . 
The  following  matrix  manipulations  can  be  found  in  [20] . 


Suppose 


(4.6) 


It  follows  that 


(4.7) 


taking  the  trace  of  both  sides 


tr  {[/>fs.[PfU  tr 


(4.8) 


^ggjWJUWUlWKX-V..^^  V,  V.VV  y.  w  ^  ^Tgygw^VC^T  77  V.’  y*  -.■  ^  ^  u_-  ^  v»  ■— 


using  aquation  (2.37)  the  terms  inside  the  trace  may  be 
rearranged 


tr  {StP‘}*  tr  {5,.,?*} 


(4.9) 


now,  multiply  both  sides  by  m* 


tr  (StP‘}A/‘*  tr  {S,_,/>‘}A/‘ 


(4.10) 


next,  a  sum  over  all  i  does  not  change  the  inequality 


i£  tr  (S.f'lw'iil  tr  {S,.,?')*/' 

•*  I- 1  *■  »- 1 


(4.11) 


similarly  add  R  to  each  side 


«• 


R  +  hZ  tr  (5t/>‘)A/‘^A  +  ^£  tr  (4.12) 


taking  the  inverse  of  both  sides  reverses  the  inequality 


U  +  jJj  tr  +  tr  (4.13) 


now  multiply  from  the  left  by  B  and  from  the  right  by  ®r 


+  tr  +  tr  {S.-.P'lA/'J  BT 


(4.14) 


i 


inequality  (4.6)  implies 


using  this  fact 


t 

l 

s;‘  +  sU*yj  tr  {S,/"|V,J  S'S 
ss;!,*srft  +  i£  tr  (s,.,/l,)wl]  bt  (4.16) 


again  taking  the  inverse  reverses  the  inequality 


tr  Brl  i 

tr  (S,.,B‘)A/'J  Brl 


(4.17) 


multiply  from  the  left  by  AT  and  from  the  right  by  A 


tr  5rJ  A * 

tr  BT 1  A 


(4.18) 


The  other  summation  term  can  be  developed  similar  to  in¬ 
equalities  (4.6)  through  (4.11),  therefore  adding  0  and 
the  other  summation  to  both  sides  does  not  affect  the 
inequality 
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C*ylr|  s;‘-fi^K  +  tr  BrJ  /l  +  jl  tr  (S,B,)U''* 

tr  (S,.l/>‘)l/‘]  Brj  tr  (S,.1f>‘)l/' 


(4.19) 


Comparing  inequality  (4.19)  with  equation  (4.5)  one  can 
see  that  the  left  side  is  just  5*.,  and  the  right  side  is 
just  5t.  Therefore  inequality  (4.6)  implies  the  following 


S  t*i  ^  S  i 


(4.20) 


Therefore,  by  inductive  reasoning,  starting  from  5.=0,  5k 
is  not  only  positive  definite,  but  it  is  also  monotonical- 
ly  increasing. 


All  that  is  left  to  do  is  to  determine  a  condition 
under  which  tho  positive  definite  monotonically  increasing 
matrix  sequence  5t  remains  bounded.  First,  assume  5k  gets 
large,  but  that  it  has  not  yet  exceeded  the  limit  given  by 
inequality  (4.21)  for  a  positive  scaler  a. 


St*  al 


(4.21) 


We  need  to  find  a  condition  such  that,  for  large  5t  , 


inequality  (4.21)  implies 


S  *♦!  ^  oi  I 


(4.22) 
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so  that  tha  limit  given  by  si  is  not  exceeded  by  subse¬ 
quent  terms  in  the  matrix  sequence.  Then,  through  induc¬ 
tive  reasoning  we  will  be  able  to  conclude  that  this 
threshold  will  never  be  exceeded.  Similar  to  the  proce¬ 
dure  followed  with  inequality  (4.6),  inequality  (4.21)  can 
be  transformed  into  inequality  (4.23) 


Q  +  +  +  ^  tr  A+l- £ 

S<2+y*rjV‘  I  +  tr  A  +  X- 


tr 

a  £  tr  (/,V‘ 

(•i 

(4.23) 


The  left  side  of  inequality  (4.23)  is  just  5*.,  .  Therefore 
if  we  set  the  right  side  of  inequality  (4.23)  less  than  or 
equal  to  al,  then  inequality  (4.22)  is  satisfied  and  s, 
will  remain  bounded  and  reach  a  steady-state.  Therefore, 
the  system  statistics  must  be  such  that  the  following 
inequality  is  satisfied 


i<x]T  tr  {r‘) l/‘$a  I 

2  t-i 


(4.24) 


Factoring  a  out  of  the  inverses 
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Si 


I 


(4.25) 


lrjVl  I  tr  M*/‘j  A 

+Q  +  ;UZ  tr  (P‘)l/‘£a  I 

2  (.1 


and 


(3  +  ay|rj^  I  +5^a'IA  +  ^|i  tr  (/>‘}A/‘j  firJ  A*^af{  tr 


(4.26) 


Dividing  both  sides  by  alpha 


a'lQ  +  AT 


I  +  fi£a'l/i  +  |£  tr  {P'}A/‘j  firJ  /,  +  ^Z  tr  {ptWl 


Z  I 

(4.27) 

As  a  was  an  arbitrarily  large  number,  « ''q  and  «-•/»  can  be 
considered  negligible,  thus  completely  eliminating  a  from 
the  equation.  If  we  also  assume 


l-t  tr  (?‘)A/‘>  0 


we  obtain 


At 


♦b[jI  tr  B'\  yl*5,I  tr  M1*" 


<  I 


(4.28) 
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This  simply  means  that  the  magnitude  of  the  eigenval¬ 
ues  of  the  left  hand  side  of  equation  (4.28)  must  be  less 
than  one.  This  is  only  a  sufficient  condition  for  the 
existence  of  a  steady-state  solution  to  St  .  As  was  stated 
before,  if  the  magnitude  of  the  maximum  eigenvalue  of  the 
left  hand  side  of  equation  (4.28)  exceeds  one,  a 
steady-state  solution  for  S*  may  still  exist,  but,  as  long 
as  the  maximum  is  less  than  one,  a  steady-state  solution 
is  guaranteed  to  exist. 


Even  though  equation  (4.28)  is  only  a  sufficient  con¬ 
dition  for  the  multivariable  case,  it  reduces  to  the  nec¬ 
essary  and  sufficient  condition  given  in  equation  (3.16) 
for  the  scalar  case.  If  we  evaluate  equation  (4.28)  for 
the  scalar  case,  then  n*l,  which  implies  n'=l.  Therefore, 
we  only  have  a  single  term  for  each  of  the  summations, 
and,  if  we  let  PJtm  U  and  PV=V  ,  then  equation  (4.28)  re¬ 
duces  to 


I  +  T  2 


(4.29) 


clearing  some  of  the  fractions 


W  a2U 


2  2  b2  +  M 


(4.30) 


adding  and  subtracting  <x’ 


-47- 


ft 


W5 


rfV-r.V.V'-W  rfVWW*..  ✓***■"  wV«r 


.r.vv-  / 


r  2) 


2  b2+M 


VT.W  *  J 


finally  this  simplifies  to  the  scalar  condition,  equation 
(3.16),  less  the  cross  term  N/ 2 


I  2)  b2  +  Z 


(4.32) 


Therefore,  as  was  stated  above,  even  though  equation 
(4.28)  is  only  a  sufficient  condition,  in  the  scalar  case, 
it  does  indeed  reduce  to  a  necessary  and  sufficient  condi¬ 
tion.  It  should  also  be  noted  that,  as  in  the  scalar 
case,  we  must  require  /»(**. ut. a,)  to  not  contain  a  purely 
additive  noise  term,  making  f*  and  ak  both  equal  to  zero. 


4 . 2  Simulationi 


In  the  following  simulations,  n  is  equal  to  two,  and 
S»  is  plotted  vs  k  for  fifty  values  of  k.  With  n  equal  to 
two,  n'  is  equal  to  three.  Therefore  if  we  let 


p-.I'  *■-!« 


then  equations  (4.2)  and  (4.3)  reduce  to 


s 


"  i  ' 

TTk-~m  tr  (S*)  (4.34) 

Inserting  equations  (4.33)  and  (4.34)  into  equation  (4.1) 
results  in 

+  +  y  tr  (Sjj  +  |  tr  (5t) 

(4.35) 

Applying  the  same  parameters  to  the  threshold  condition, 
equation  (4.28),  results  in 

wl  +/tr£  I  A<  I  (4.36) 

and  similar  to  the  scalar  case,  the  optimal  cost  to  go 
from  time  0  to  time  k  is 

Jk-xlSkx  o  (4.37) 

The  first  four  figures  show  the  evolution  of  S* , 
equation  (4.35),  for  different  values  of  m  and  w.  The 
maximum  eigenvalue  of  the  threshold  condition  is  shown  at 
the  top  of  each  figure. 


‘  rj»  ”  »  “  ■  ^  Mi  wi  wit  w  «  \r»  i/v  L/v\N  '  \  •  v  ».v.  *  j  fj  VKKTvk  rjr  ru*  YJ rW irjTrjXVwVWYFV  KVTYTW 


Also,  with  *  equal  to  two,  s,  is  a  two  by  two  matrix, 
and,  as  it  is  symmetric,  5»(1,2)  *  st(2,l).  Therefore  only 
three  plots  appear  in  each  figure.  The  next  four  figures 
show  the  respective  evolution  of  the  optimal  cost  given  in 
equation  (4.37),  with  x»  equal  to  a  unit  vector. 


i 


Unlike  the  scalar  case,  we  do  not  have  large  varia¬ 
tions  in  magnitude  within  each  individual  figure  so  there 
is  no  reason  to  take  the  logarithm.  Therefore,  the  actual 
values  of  S,  and  J»  are  plotted  vs  k,  but  one  should  note 
that  there  are  still  large  variations  in  magnitude  between 
the  figures.  The  following  values  apply  to  all  simula¬ 
tions.  Jl=0 . 1 ,  and 


I 

a 

o 

0.81 

c-f0 

o.ol 

sJ0-1! 

r 

Lo.o 

0.2J 

v  lo.o 

1  .oj 

L0.2J 

Sk  (1,2) 


EVOLUTION  OF  Jk 

LAMBDA  MAX  -  1.163 


Q 


r-: 


ij 


CHAPTER  5 


V 


CONCLUSIONS 

We  began  by  introducing  a  very  general  system  de¬ 
scription.  We  showed  that  this  system  description  can 
represent  anything  from  a  simple  deterministic  system  to  a 
very  complex  nonlinear  stochastic  system,  and  gave  several 
examples.  After  giving  the  optimal  finite-horizon  solu¬ 
tion  for  several  known  simpler  cases,  we  then  stated  the 
finite-horizon  optimal  controller  for  the  proposed  multi- 
variable  nonlinear  stochastic  system,  subsequently  proving 
the  general  result . 

The  practicality  of  the  infinite-horizon  control  ap¬ 
proach  was  pointed  out,  prompting  an  investigation  into 
the  steady-state  characteristics  of  the  Riccati-like  equa¬ 
tions  involved  in  the  finite-horizon  control  solutions.  A 
necessary  and  sufficient  condition  for  the  existence  of  a 
steady-state  solution  was  established  for  the  scalar  case. 
Due  to  the  complexities  involved  in  the  matrix  case,  only 
a  sufficient  condition  was  developed  that  guaranteed  the 
existence  of  a  steady-state  solution.  This  sufficient 
condition  for  the  multivariable  case  was  shown  to  reduce 
to  the  necessary  and  sufficient  condition  for  the  scalar 
case.  Extensive  simulations  were  provided  for  both  cases 
to  verify  these  conditions. 
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Possible  future  work  in  this  area  might  investigate 
the  continuous-time  case.  The  continuous-time  counterpart 
of  the  general  controller  proposed  in  this  work  would  most 
likely  be  needlessly  complex.  Therefore,  it  would  seem 
advantageous  to  develop  a  sampling  procedure  for  continu¬ 
ous-time  systems,  such  that,  the  resulting  discretized 
system  would  fit  into  the  framework  of  the  discrete-time 
solution  already  given  in  this  work. 
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************************* 

*  MARK  CIANCETTA  * 

*  20  JUL  87  * 

*  THIS  PROGRAM  COMPUTES  THE  THRESHOLD  VALUE,  * 

*  AND  SK  FOR  THE  SCALAR  CASE  * 

************************* 

REAL  M,N 


A  *  1.1 
B  *  1.0 
Q  =  1.0 
T  =  .2 
SK  *  0. 

R  *  1.0 

*  *  *  *  INPUT  VALUES  FOR  M,  N,  AND  W  *  *  *  * 

WRITE  (5,5) 

FORMAT  (IX, 'M  =  ?') 

READ  (5,10)  M 
FORMAT  (F5.2) 

WRITE  (5,15) 

FORMAT  (IX, 'W  =  ?') 

READ  (5,10)  W 
WRITE  (5,17) 

FORMAT  (IX, 'N  *  ?') 

READ  (5,10)  N 

*  *  *  *  CALCULATE  THE  THRESHOLD  CONDITION  *  *  *  * 
g  «  (A* *2  +  W/2)  -  ( (A*B  +  n/2) **2/ (B**2  +  M/2)) 

*  *  *  *  OUTPUT  VALUE  AND  PROMPT  TO  CONTINUE  *  *  *  * 
WRITE  (5 ,20)  g 

FORMAT  (IX, 'THE  THRESHOLD  VALUE  g  ■  ',F6.3) 

WRITE  (5,25) 

FORMAT  (IX, 'ENTER  1  TO  CONTINUE') 

READ  (5,30)  I 
FORMAT  (12) 

IF  (I.NE.l)  GO  TO  50 

*  *  *  *  0PEN  A  FILE  AnD  INITIALIZE  IT  *  *  *  * 

OPEN  (UNIT*10,FILE=' SK' , STATUS*' NEW' ) 

WRITE  (5,35)  g 
WRITE  (10,35)  g 
WRITE  (5,35)  SK 
WRITE  (10,35)  SK 
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************************** 

*  MARK  CIANCETTA  * 

*  19  JUN  87  * 

*  THIS  PROGRAM  COMPUTES  THE  THRESHOLD  VALUE,  SK,  * 

*  AND  JK  FOR  THE  MULTIVARIABLE  CASE  * 

************************** 

DIMENSION  A  (2 ,2)  ,AT(2,2)  ,B(2)  , BBT  (2, 2)  ,Z(2,2) 
DIMENSION  Q(2 ,2)  ,SKT(2,2)  ,SK(50,4)  ,ZI(2,2) 

REAL  M,  JK(50) 

DATA  A/. 1,0., .8, . 2/ , B/ . 1 , .2/ ,SK/200*0./ , 

*  Q/l.,0.,0.,1./, JK/ 50*0./ 

R  =  .1 

*  *  *  *  INPUT  in  AND  w  *  *  *  * 

WRITE  (5,5) 

FORMAT  (IX, 'm  =  ?') 

READ  (5,10)  M 
FORMAT  (F5.2) 

WRITE  (5,15) 

FORMAT  (IX, 'w  =  ?') 

READ  (5,10)  W 

*  *  *  *  CALCULATE  AT(ranspose)  *  *  *  * 

AT  (1,1)  =  A(1 , 1) 

AT (1 , 2)  *  A (2 , 1) 

AT  (2 , 1)  *  A(1 ,2) 

AT  (2, 2)  *  A  (2 ,2) 

*  *  *  *  CALCULATE  BBT (ranspoie)  *  *  *  * 

BBT  (1,1)  =  B(l)  **2 
BBT (1,2)  =  B(l) *B (2) 

BBT  (2,1)  =  B(l)  *B(2) 

BBT  (2, 2)  =  B  (2)  **2 

*  *  *  *  CALCULATE  Z  =  I  +  1/m  BBT  *  *  *  * 

Z(l,l)  =  1  +  BBT  (1,1)  /M 
Z  (1 , 2)  =  BBT (1 , 2) /M 
Z  (2 , 1)  =  BBT  (2 , 1)  /M 
Z  (2 , 2)  =  1  +  BBT (2, 2) /M 

*  *  *  *  det .  OF  Z  *  *  *  * 

D  =  Z  (1 , 1)  *Z  (2 , 2)  -  Z(1 ,2)  *Z(2,1) 

*  *  *  *  CALCULATE  INVERSE  *  *  *  * 
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ZI(l.l) 

zia,2) 

ZI(2,1) 
ZI (2 ,2) 


Z  <2 ,2)  /D 
-Z(l,2)/D 
-Z(2,l)/D 
Z(l,l)/D 


*  *  *  *  CALCULATE  AT (ranspose) ZI (nverse) A  ->  Z  *  *  * 


CALL  MULT  (AT ,  ZI ,  Z) 

CALL  MULT (Z ,A, Z) 

*  *  *  *  ADD  wl  *  *  *  * 

Z(l,l)  =  Z(l,l)  +  w 
Z  (2 , 2)  »  Z  (2 ,2)  +  W 

*  *  *  *  CALCULATE  MAG  OF  THE  EIGENVALUES  *  *  *  * 

r  -  7  n  n  +  7  O  9) 

D  =  4*  (Z  (1 , 1)  *Z  (2 , 2)  -  Z  (1 ,2)  *Z  (2 , 1)  ) 

D  =  C**2  -  D 

IF  (D . GE . 0 . ) THEN 

El  =  ABS ( (C  +  SORT (D) )  /  2) 

E2  =  ABS ( (C  -  SORT  (D) ) /2) 

ELSE 

D  ■  -D 

El  *  SORT (C**2  +  D)/2 
E2  =  SORT  (C* *2  ~  D>/2 
ENDIF 

*  *  *  *  FIND  LARGEST  EIGENVALUE  *  *  *  * 

IF  (El .GT.E2)THEN 
E  =  El 
ELSE 

E  =  E2 
ENDIF 

*  *  *  *  OUTPUT  THRESHOLD  VALUE  AND  CONTINUE?  *  *  *  * 
WRITE  (5,30)  E 

FORMAT  (IX, 'THE  MAXIMUM  EIGENVALUE  IS  ',F6.3) 

WRITE  (5,32) 

FORMAT  (IX, 'ENTER  1  TO  CONTINUE’) 

READ  (5,33)  N 
FORMAT  (12) 

IF  (N.NE.l)  GO  TO  50 

*  *  *  *  OPEN  FILE  AND  OUTPUT  S(0)  AND  J(0)  *  *  *  * 
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oo  non  non  noo  noo  non  oon  ooo  ooo 


OPEN  (UNIT=10,FILE=' SK' , STATUS-' NEW' ) 

WRITE  (5,35)  SK(1,1)  ,SK(1,2)  ,SK(1,3)  ,SK<1,4)  ,JK(1) 
WRITE  (10,35)  SK(1,1)  ,SK(1,2)  ,SK(1,3)  ,SK(1,4)  ,JK(1) 

*  *  *  *  STORE  PREVIOUS  SK  IN  SKT(emp)  *  *  *  * 

DO  40  1=1,50 

SKT (1,1)  =  SK (I , 1) 

SKT (1 ,2)  =  SK (1 , 2) 

SKT  (2,1)  =  SK  (1 , 3) 

SKT (2 ,2)  =  SK(I ,4) 

*  *  *  *  CALCULATE  tr  [SKT]  /2  *  *  *  * 

TR  =  (SKT (1,1)  +  SKT (2 ,2) ) /2 

*  *  *  *  CALCULATE  B (transpose) SKB  *  *  *  * 

BTSKB  =  B(l)  **2*SKT  (1,1)  +  B(l)  *B(2>  *SKT  (2 , 1) 
*  +  B(1)*B(2)  *SKT  (1,2)  +  B(2)  **2*SKT  (2 ,2) 

*  *  *  *  THE  INVERSE  ASSOCIATED  WITH  S(k+1)  *  *  *  * 

C  =  1/(R  +  BTSKB  +  M*TR) 

CALL  MULT (SKT, BBT, Z) 

CALL  MULT (Z , SKT , Z) 

*  *  *  ♦  SUBTRACT  THIS  FROM  SK  *  *  *  * 

SKT  (1,1)  =  SKT (1,1)  -  C*Z(1,1) 

SKT (1 , 2)  =  SKT (1 , 2)  -  C*Z(1,2) 

SKT (2,1)  =  SKT (2,1)  -  C*Z(2,1) 

SKT (2 ,2)  =  SKT (2 ,2)  -  C*Z(2,2) 

*  *  *  *  MULT  FROM  THE  LEFT  BY  A (transpose)  *  *  *  * 

CALL  MULT (AT , SKT , Z) 

*  *  *  *  MULT  FROM  THE  RIGHT  BY  A  *  *  *  * 


CALL  MULT (Z , A, Z) 


*  *  *  *  ADD  Q  AND 

(w/2) 

*  tr  (SK)  * 

SK (1+1 , 1)  = 

Z(l,l) 

+  1.  +  W*TR 

SK (1+1 ,4)  = 

Z(2 ,2) 

+  1 .  +  W*TR 

SK (1  +  1 ,2)  = 

Z (1 , 2) 

SK(I  +  1 , 3)  = 

Z  (2 , 1) 

*  *  *  *  CALCULATE 

JK  *  * 

*  * 
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JKCI+1)  =  (SK(I  +  1,1)+SK  (1+1,2)  + 

*  SK(I  +  l,3)+SK(I+l,4))/2 

*  *  *  *  OUTPUT  THE  RESULTS  *  *  *  * 

J  =  1+1 

WRITE (5 , 35) SK ( J , 1) ,SK(J,2) ,SK(J,3) ,SK(J,4) ,JK(J) 
WRITE (10 , 35)  SK(J ,  1)  ,SK(J,2)  ,SK(J,3)  ,SK(J,4)  ,JK(J) 
5  FORMAT  (IX ,5  (2X , F15 . 5) ) 

0  CONTINUE 

STOP 
END 
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************************ 

*  MARK  CIANCETTA  * 

*  15  JUNE  1987  * 

*  THIS  PROGRAM  MULTIPLIES  TWO  2X2  MATRICES,  * 

*  A  AND  B,  AND  RETURNS  THE  PRODUCT  IN  C  * 

************************ 

SUBROUTINE  MULT(A,B,C) 

DIMENSION  A (2, 2)  ,B(2,2)  ,C(2,2) 

DO  4  I  =  1,2 

DO  3  J  =  1,2 

C(I,J)  =  A(1 , 1)  *B(1 ,  J)  +  A(1 ,2)  *B(2  ,  J) 
CONTINUE 
CONTINUE 

RETURN 


